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We propose a simple and general procedure based on a recently introduced approach that uses an 
importance-sampling Monte Carlo algorithm in the disorder to probe to high precision the tails of 
ground-state energy distributions of disordered systems. Our approach requires an estimate of the 
ground-state energy distribution as a guiding function which can be obtained from simple-sampling 
simulations. In order to illustrate the algorithm, we compute the ground-state energy distribution 
of the Sherrington-Kirkpatrick mean-field Ising spin glass to eighteen orders of magnitude. We find 
that if the ground-state energy distribution in the thermodynamic limit is described by a modified 
Gumbel distribution as previously predicted, then the value of the slope parameter m is clearly 
larger than 6 and of the order 11. 

PACS numbers: 75.50.Lk, 75.40.Mg, 05.50.+q 



I. INTRODUCTION 

The ground-state energy distributions of disordered 
systems and in particular of spin glasses 0, 0, 0, 01 have 
recently received considerable attention from different 
groups. Bouchaud et al. [j| have studied the ground- 
state energy distribution of short-range spin glasses in 
two and three dimensions and find them to become Gaus- 
sian in the thermodynamic limit 6], while other groups 
III El nn d that the ground-state energy distri- 
bution of the mean-field Sherrington-Kirkpatrick (SK) 
model 1,1] remains skewed in the thermodynamic limit 
and apparently can be described by a modified Gumbel 
distribution In a recent paper, we have been able 

to show that this transition from a Gaussian to a non- 
Gaussian limiting distribution coincides with the tran- 
sition from long-range to infinite-range behavior for a 
one-dimensional spin glass with long-range power law in- 
teractions For uncorrelated random variables, when 
the power of the exponent in the aforementioned model is 
large, the central limit theorem holds and a Gaussian dis- 
tribution can be expected. But when the interactions are 
infinite-ranged and the variables are strongly correlated, 
such as in the case of the SK model, deviations from the 
central limit theorem can be expected, and, as previously 
mentioned, the data seem to follow a modified Gumbel 
distribution with a slope parameter expected to be 
m ~ 6 [see Eq. (|10f) ]. This result is rather puzzling: What 
is the relationship between the ground-state energy dis- 
tribution of the SK model and extreme value statistics? 
Recently Bertin and Clusel [14] analytically derived the 
relationship between extreme value statistics and random 
sums of correlated variables. Thus it might be plausible 
that similar arguments could be applied to the SK model 
in order to explain the skewed energy distributions. In 
particular, the results of Bertin and Clusel show for sums 



of correlated variables that there is not necessarily an un- 
derlying extreme process in order to obtain a modified ex- 
treme value distribution. This means that the parameter 
m, expected to be integer when studied from the context 
of extreme value statistics, can also be noninteger, as has 
been found by by Bramwell et al. [TH fl(| when studying 
large-scale critical fluctuations in correlated systems. 

Traditionally El, simple-sampling techniques, 
where iV samp disorder realizations are computed and sub- 
sequently binned in order to obtain the ground-state 
energy distribution, are used. If iV samp samples are 
computed, then the maximal "resolution" of a bin is 
~ 1/-/V samp , and thus ~ 10 6 samples have to be com- 
puted to resolve six orders of magnitude in the histogram. 
Therefore, the probing of tails in ground-state energy 
distributions becomes quickly intractable making it dif- 
ficult to determine if a given fitting function properly 
describes the tails of a ground-state energy distribution, 
as is the case for the SK model. Multicanonical methods 
[TtIITsI ] have been used before to overcome the limitations 
of simple-sampling techniques in order to probe tails of 
overlap distribution functions in spin glasses 0,|2(j. In 
this paper we outline a simple algorithm related to multi- 
canonical approaches based on ideas presented in Ref. |^] 
that also overcomes the limitations of simple-sampling 
techniques by performing an importance-sampling sim- 
ulation of the ground-state energy distribution in the 
disorder with a guiding function computed from simple- 
sampling simulations. Similar approaches have been used 
before in the studies of distributions of sequence align- 
ment scores pi) , free-energy barriers in the Sherrington- 
Kirkpatrick model , as well as fluctuations in classical 
magnets [23| (albeit the latter without disorder). 

By computing the tails of the disorder distribution of 
the SK model to up to 18 orders of magnitude we show 
that a modified Gumbel distribution fits the data well, 
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yet with small systematic deviations. If the ground-state 
energy distribution in the thermodynamic limit is a modi- 
fied Gumbel distribution, then the slope parameter of the 
Gumbcl distribution is considerably larger (m « 11) than 
found in previous studies using simple-sampling tech- 
niques 0- 

We outline the method in Section [D] and apply it to 
the ground-state energy distribution of the SK model in 
Section IIIII We conclude in Sec. II VI and discuss the im- 
plications of our results on studies of ground-state energy 
distributions. 



II. SIMULATION OF GROUND-STATE 
ENERGY DISTRIBUTIONS 

A disordered system is defined by a Hamiltonian 
7ij(C), where the disorder configuration J is chosen 
from a probability distribution P(J) and C denotes the 
phase-space configuration of the system. The ground- 
state energy E of a given disorder configuration J is 
defined by 

E(J) = mm H j (C). (1) 

Together with the disorder distribution V(J), this de- 
fines the ground-state energy distribution 

P{E) = J dJV(J)S[E-E{J)}. (2) 

To study the thermodynamic limit of the ground-state 
energy distribution, we use the standardized form P s (x) 
defined by the equation 

P (E) = ± Ps ( , ( 3 ) 

where [£] av and oe = {[E 2 ] &v — [E]^) 1 ^ 2 are the average 
and the standard deviation of the ground-state energy, 
respectively [24|. Here [• • -] av represents an average over 
the disorder. 

A. Simple sampling 

Because P(E) cannot be determined analytically for 
most disordered systems, approximative methods such 
as Monte Carlo simulations have to be used. A standard 
approach is to use a simple-sampling Monte Carlo algo- 
rithm to study P{E). iVgamp independent disorder config- 
urations Ji are chosen from P(J) and the ground-state 
energy is calculated for each disorder configuration. The 
calculation of the ground-state energy in itself is a diffi- 
cult optimization problem and can often only be solved 
approximatively |25, 26] . From the ground-state energies 
of these disorder configurations, the ground-state energy 
distribution can be estimated as 

JV samp 

P(E) = —- ]T 6[E-E(J i )}, (4) 



so that the averages of functions with respect to the dis- 
order are replaced by averages with respect to the N sainp 
random samples. The functional form of the ground-state 
energy distribution and its parameters can for example 
be estimated by a maximum likelihood fit of a distribu- 
tion Fg(E) with parameters 9 to the data |23|. It cor- 
responds to the distribution for which the observed data 
have the largest probability. Due to the limited range 
of energies sampled by the simple-sampling algorithm it 
is often difficult or even impossible to quantify how well 
the tails of the distribution are described by a maximum- 
likelihood fit, and improving the data by increasing the 
number of random samples is generally computationally 
very expensive. Therefore other methods have to be used 
to study the tails of the distribution, and we outline in 
the next section a simple method that can be applied if 
a good analytic fit to the data can be found. 



B. Importance sampling with a guiding function 

Assuming that we find a function Fg(E) which ac- 
curately describes the ground-state energy distribution 
as calculated with a simple-sampling simulation, we now 
proceed by sampling the ground-state distribution with 
an importance-sampling Monte Carlo algorithm in the 
disorder 0, 12^, Hij and use Fg (E) as a guiding func- 
tion |3(il ] . We start from a random disorder configuration 
J = Jo with ground-state energy E(J ). From the z-th 
configuration J7i, we generate the i + 1-th configuration 
Ji+i by the following Metropolis-type [3l| update: 

1. Choose a candidate disorder configuration J' by 
replacing a subset of J chosen at random (e.g., a 
single bond chosen at random) with values chosen 
according to P(J) [this requires that P(J) can be 
written in a product form] and calculate its ground- 
state energy E(J'). 

2. Set Ji+i = J' with probability 

p . f Fg [E(Jj)] \ 

Accept = mm I ^rj^jjy^ > j (5) 

and Ji+\ = Ji otherwise. 

Using the importance-sampling Monte Carlo algorithm, 
a disorder configuration J is visited with probabil- 
ity l/Fg[E(J)), such that the probability to visit a 
disorder configuration with ground-state energy E is 
P(E)/Fg(E). If F e (E) = P(E), then each energy is 
visited with the same probability resulting in a flat- 
histogram sampling of the ground-state energy distribu- 
tion. To be able to study a finite range of energies and 
to avoid the trapping of the algorithm in an extremal re- 
gion of the energy space, the range of energies that the 
algorithm is allowed to visit can be restricted. 

The main difference to the simple-sampling algorithm 
described in Section [II Al is that successive configurations 
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visited by the algorithm are not independent. There- 
fore an analysis of the results requires a quantification 
of the correlations among configurations visited by the 
algorithm. This is a standard problem in Markov Chain 
Monte Carlo simulations, and can for example be solved 
with the help of the exponential autocorrelation time r 
of the energy, which is the number of Monte Carlo steps 
after which the autocorrelation function of the ground- 
state energy 



( E i E i+Ai) ~ ( E i) ( E i+Ai) 



(0) 



decays to 1/e |2g|. Here Ei is the ground state energy 
after the i-th. Monte Carlo step and (. . .} refers to the av- 
erage over Monte Carlo time. In order to ensure that the 
visited ground-state configurations are not correlated, we 
only use every 4r-th measurement. Once the autocorre- 
lation effects have been quantified, the data can be an- 
alyzed with the same methods as the simple-sampling 
results [see Eq. @] ■ 



III. APPLICATION: 
SHERRINGTON-KIRKPATRICK MODEL 



TABLE I: Maximum-likelihood fit to the simple-samp ling 
Monte Carlo data calculated for the SK model in Ref. Ha . 
For each system size N, iV sam p = 10 J random samples have 
been generated. The maximum-likelihood fit has been per- 
formed by binning the data into 50 bins. The error bars and 
the parametric estimates of [£]„ and o\e [see Eqs. 1 1 1 li and 
11211 1 have been generated by a bootstrap method (see text). 
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[Slav" 


_ a 
<JE 


16 


-10.373(6) 


2.48(4) 


4.92(16) 


-10.634(4) 


1.180(3) 


24 


-16.189(6) 


2.81(5) 


4.99(17) 


-16.481(4) 


1.326(3) 


32 


-22.076(6) 


3.23(6) 


5.59(19) 


-22.375(5) 


1.432(3) 


48 


-33.935(8) 


3.81(8) 


6.24(26) 


-34.249(5) 


1.590(4) 


64 


-45.848(10) 


3.97(9) 


5.87(26) 


-46.196(5) 


1.712(4) 


96 


-69.840(9) 


4.53(9) 


6.37(25) 


-70.205(6) 


1.867(5) 


128 


-93.918(13) 


4.89(14) 


6.49(39) 


-94.305(6) 


1.997(5) 



"Parametric estimate calculated from fi, u, and m. 



studying the I0 5 ground-state energies calculated for dif- 
ferent system sizes in Ref. We bin the data into 50 
bins and perform a maximum-likelihood fit with a mod- 
ified Gumbel distribution pi LL2J which seems to fit the 
simple-sampling data well 



The Sherrington-Kirkpatrick 1 1 11] model is defined by 
the Hamiltonian 



(7) 



i<3 



where the S% = ±1 (i = l,...,N) are Ising spins, 
and the bonds J = {Jij} are identically and indepen- 
dently distributed random variables chosen from a Nor- 
mal distribution with zero mean and standard deviation 
(N — l) -1 / 2 . The sum is over all spins in the system. 

The ground-state energy of a given disorder configura- 
tion J is defined as 



E(J) = min Hj({Si}). 



(8) 



We are interested in the thermodynamic limit of the 
ground-state energy distribution and define the limiting 
distribution Pqo(x) in analogy to Eq. © by 



P{E) . _Lp„ (£z!S. 

OE \ &E 



(9) 



For the SK model several optimization algorithms, such 
es extremal optimization |32j |. hysteretic optimization 
[3^| . as well as other algorithms such as genetic and 
Bay esian al gori thms |25ll26l |. and even parallel tempering 
1 1.1 I34I I'iol l36| are available. In this work we compute 
the ground states of the system for the sake of simplicity 
using parallel tempering Monte Carlo |3]|. Note that a 
combination of the proposed importance-sampling sim- 
ulation in the disorder with other more efficient algo- 
rithms might yield more detailed results. We start by 



G^,v, m {E) K exp 



E-n [E-u 
m m exp 



(10) 



The modified Gumbel distribution is parameterized by 
the "location" parameter fi, the "width" parameter 1/, 
and the "slope" parameter m. To obtain error bars we 
generate 200 bootstrap replicates j^] of the data and 
perform the maximum-likelihood fit for each bootstrap 
replicate. Because the average and standard deviation 
of the modified Gumbel distribution are related to the 
parameters /i, v, and m by 



[E]ay = fl + E m V 



(11) 
(12) 



where E m and <j m are the average and standard devi- 
ation of Go,i, m , a parametric estimate of the average 
energy and its standard deviation can be calculated for 
the data (and each bootstrap replicate). The results of 
the analysis are summarized in Table [I] Note that for 
the simple-sampling data we find to « 6 as suggested in 
Ref. |3 and that the parametric estimates of [E] av and ue 
agree within error bars with the direct estimates calcu- 
lated in Ref. using Eqs. (4) and (5) therein. 



Importance-sampling simulation 

Using the maximum-likelihood estimates shown in Ta- 
ble[IJ we now perform an importance-sampling simulation 
in the disorder as described in Section III Bl To perform 
a step in the Monte Carlo algorithm, we choose a site 
at random, replace all bonds connected to this site (the 
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TABLE II: Input parameters for the guiding function used in 
the simulation. Initial runs using the parameters shown in 
Table H] indicated that m > 6, so that we have chosen m = 8 
for the production runs. The parameters fj, and v have been 
calculated from [-E] av , o e , and m with the help of Eqs. 1111 
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FIG. 1: Autocorrelation function as defined in Eq. JSJ for 
system sizes N — 16 (circles) and N — 128 (triangles) for the 
simulation with parameters shown in Table IH1 The value 1/e 
is marked by the horizontal dotted line (Ai is measured in 
Monte Carlo steps). 



expected change in the ground-state energy is then of 
the order ~ 1/-/V), calculate the ground-state energy of 
the new configuration, and accept the new configuration 
with the probability given in Eq. JS}, which in this case 
is a modified Gumbel distribution, Eq. (|10|) . with the pa- 
rameters listed in Table |n] To avoid a trapping of the 
simulation in the double-exponential forward tail of the 
distribution, we limit the maximum energy allowed in the 
simulation by [E] av + 3cte. Initial runs indicated that the 
estimates of m shown in Table [I] are too small, and we 
therefore choose an estimate of m = 8 for the production 
runs, /i and v are determined from the simple-sampling 
results for the average and standard deviation with the 
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FIG. 2: Unsealed ground-state energy distributions of the 
Sherrington-Kirkpatrick model for different system sizes, ob- 
tained by the guiding-function simulation with the parameters 
given in Table ITT1 For each system size, between 92 and 686 
independent samples are simulated. 



TABLE III: Three-parameter fit in the parameters /i, v, and 
m to the data [rescaled to x — (E — [E]a. v ) /cte] for the SK 
model. For each system size N, between JV samp = 92 and 686 
independent samples have been generated, z — vjm describes 
the asymptotic behavior of the single-exponential tail (error 
bar obtained from independent fits). y 2 /dof represents the 
X 2 per degree of freedom of the fit |24j|. 
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help of Eqs. 1)11(1 and l)12[i . and a summary of the input 
parameters used is shown in Table ITT1 It is important to 
note that this change in the parameters does not lead to a 
systematic error or bias in the results and merely consti- 
tutes a change of the guiding function, which can either 
improve or degrade the range of energies visited by the 
algorithm. Figure ^ shows the energy-energy autocorre- 
lation function for system sizes N = 16 and 128. Auto- 
correlation times are of the order of 400 to 700 Monte 
Carlo steps resulting in 92 to 686 independent measure- 
ments for the different system sizes. While the number 
of samples used is small, the method is able to probe the 
tails in this particular case down to 18 orders of magni- 
tude, a result impossible to obtain with simple-sampling 
techniques. Figure|21shows the unsealed ground-state en- 
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FIG. 3: Scaled ground-state energy distributions of the 
Sherrington-Kirkpatrick model for different system sizes. The 
curves are scaled by the average energy and standard devia- 
tion obtained from simple-sampling simulations shown in Ta- 
ble U] The scaled distributions obtained for different system 
sizes collapse on a common curve, although finite size effects 
are clearly visible in the tail. The inset shows the behavior of 
the ground-state energy distribution at x — — 11.5 (vertical 
dotted line in the main panel) as a function of the system size 
iV in order to illustrate finite-size effects. 



ergy distributions obtained from the simulation. Figure|3| 
shows the same distributions standardized with the aver- 
age energy and the standard deviation as measured from 
simple-sampling simulations. We choose to standardize 
with [E] av and erg from the simple-sampling simulations 
because of the greater number of statistically - for these 
quantities - relevant measurements. The standardized 
distributions obtained for different system sizes collapse 
on a single curve, although finite size effects are clearly 
visible in the tail of the distribution (see also inset of 
Figure OJl • 

Next we try to determine the functional form of the dis- 
tribution by fitting a modified Gumbcl distribution to the 
data. We have used the standard fitting function of the 
gnuplot package, which implements the nonlinear least- 
squares (NLLS) Marquardt-Levenberg algorithm. The 
results of the fit are shown as /z, v, and m in Table ITTT1 
In particular, the values for m are approximately two 
times larger than previous results Q- 

Figure0]shows the ground-state energy distribution for 
N = 128 together with the resulting fit. The data have 
been standardized by the simple-sampling results and fit- 
ted with a modified Gumbel distribution. To examine the 
systematic deviations more closely, the inset shows the 
deviation eg t of the fit from the observed data scaled by 
the measurement error AP obtained from a bootstrap 
analysis. If no systematic deviations were present, one 
would expect the deviations to be of order one with an 



FIG. 4: Scaled ground-state energy distribution for N = 128 
and the fit to the modified Gumbel distribution as described 
in the text. The inset shows ent = (Put — f 128 ) / AP128 , the fit 
deviation normalized by the error of the calculated ground- 
state energy distribution. 



irregularly changing sign. Indeed eat varies between — 1 
and 1, showing the good quality of the fit. Nevertheless, 
the differences between fit and data have clear systematic 
deviations, which might suggest that a modified Gum- 
bel distribution is not the correct asymptotic probability 
distribution function. Note that we also excluded by a 
visual comparison that the distribution has the form of 
a Tracy- Widom distribution [H, |3!j • 

Because the systematic deviations between fitted func- 
tion and data decrease with system size (not shown) 
and the fitted function agrees relatively well with the 
data over several orders of magnitude, it can be sur- 
mised empirically that a modified Gumbel distribution 
might present a good description of the limiting distri- 
bution function for the SK model. Thus we attempt 
to determine the asymptotic behavior of the tail by ex- 
trapolating to infinite system size by fitting a power law 
m(N) = Woo + aN~ b to the data for the values of m 
shown in Table ITTT1 The resulting limiting value is 

m M = 10.9(5), (13) 

where for the fit % 2 /dof = 0.331. The data points to- 
gether with the resulting fit are shown in Fig. [S] The in- 
set shows a similar analysis for z = v/m which describes 
the behavior in the exponential tail of the Gumbel dis- 
tribution, resulting in Zoa = 0.34(4). 



IV. CONCLUSIONS 

In this paper we have explained an importance- 
sampling algorithm with a guiding function to simulate 
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N 

FIG. 5: Finite-size dependence of the parameter m [see 
Eq. llUIl l from the fits to the modified Gumbel distribution, 
together with a power-law fit (see text). The inset shows a 
similar analysis for the tail parameter z = u/m. The data 
converge to moo = 10.9(5) and Zoo = 0.34(4), respectively. 

the ground-state energy distribution of a disordered sys- 
tem to high order using a considerably smaller numerical 
effort than with simple-sampling techniques. When com- 
pared to full multicanonical simulation schemes to sam- 
ple distributions, such as used in Ref. and Ref. |2^, 
our algorithm has several advantages due to its simplic- 
ity: Instead of iterating towards a good guiding function, 
which may be quite expensive computationally, we use a 
maximum likelihood fit as a guiding function. Therefore, 
the algorithm we propose is straightforward to imple- 
ment and considerably more efficient than traditional ap- 
proaches, provided a good guiding function, i.e., a good 
maximum-likelihood fit to the simple-sampling results, 
can be found. Note also that the method can be gen- 
eralized to any distribution function, such as an order- 
parameter distribution. 



We have illustrated the algorithm with the compu- 
tation of the ground-state energy distribution of the 
Sherrington-Kirkpatrick model and find that the ground- 
state energy distribution can be described over several 
orders of magnitude by a modified Gumbel distribution 
albeit with systematic deviations. Therefore, if the limit- 
ing probability distribution (N — > oo) of the SK model is 
a modified Gumbel distribution, it has a slope parameter 
to « 11, a value significantly larger than estimated be- 
fore. Simulations with more efficient ground-state search 
algorithms in order to probe larger system sizes would be 
desirable in order to see if the aforementioned systematic 
deviations become negligible for large N. Note that our 
results for the mean energy as well as the fluctuations 
o~e are not influenced by the importance sampling tech- 
nique because the method probes the tails where proba- 
bilities are small and thus contributions to the moments 
are negligible. 

Since the method can be applied more generally to 
systems where simple-sampling results exist, revisiting 
the one-dimensional Ising chain with random power- 
law interactions |ljj together with a better ground-state 
search algorithm would be desirable in order to probe the 
crossover from mean-field to non-mean-field behavior in 
more detail. 
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